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ON THE DEFINITION OF SURFACE POTENTIALS FOR FINITE-DIFFERENCE 

OPERATORS* 

S. V. TSYNKOV* 


Abstract. For a class of linear constant-coefficient finite-difference operators of the second order, we 
introduce the concepts similar to those of conventional single- and double-layer potentials for differential 
operators. The discrete potentials are defined completely independently of any notion related to the ap- 
proximation of the continuous potentials on the grid. We rather use an approach based on differentiating, 
and then inverting the differentiation of, a function with surface discontinuity of a particular kind, which 
is the most general way of introducing surface potentials in the theory of distributions. The resulting 
finite-difference “surface” potentials appear to be solutions of the corresponding system of linear algebraic 
equations driven by special source terms. The properties of the discrete potentials in many respects resemble 
those of the corresponding continuous potentials. Primarily, this pertains to the possibility of representing a 
given solution to the homogeneous equation on the domain as a variety of surface potentials w ith the density 
defined on the domain’s boundary. At the same time, the discrete surface potentials can be intei preted as one 
specific realization of the generalized potentials of Calderon’s type, and consequently, their approximation 
properties can be studied independently in the framework of the difference potentials method by Ryaben kii. 
The motivation for introducing and analyzing the discrete surface potentials was provided by the problems 
of active shielding and control of sound, in which the aforementioned source terms that drive the potentials 
are interpreted as the acoustic control sources that cancel out the uivwanted noise on a predetermined region 
of interest. 

Key words, linear differential operator, discontinuous solution, distribution, monopole, dipole, fun- 
damental solution, convolution, boundary integral, single- and double-layer potentials, surface density of 
the potential, linear difference operator, multi-layer grid boundary, discontinuous grid function, discrete 
monopole and dipole layers, Calderon’s potential, difference potential 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. For the systems of linear algebraic equations that originate from mesh discretizations 
of elliptic differential operators, we define solutions of the special structure that may be considered immediate 
discrete analogues of the continuous single- and double-layer potentials, i.e., surface potentials with the 
monopole and dipole type density, respectively. Our definition reproduces, on the mesh level, the following 
general approach to introducing surface potentials of linear constant-coefficient differential operators that 
is common in the theory of distributions. For a given domain of interest and its complement, and for a 
given linear homogeneous differential equation, the potential is defined as its piece-wise smooth classical 
solution that has, generally speaking, a discontinuity along the interface separating the domains. Applying 
the corresponding differential operator to this solution throughout the entire region of interest, we generate a 
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particular collect ion of singular layers (6-type, distributions) along the interface. In the case of second-order 
differential operators, only the discontinuity in the function itself and its first derivatives will matter, and the 
corresponding singular layers on the interface will only be of the monopole and dipole type (i.e., ,5 function 
and its first derivatives). Subsequently, the original piece-wise smooth solution on both subdomains can be 
reconstructed as convolution of the resulting singular layers with the fundamental solution of the differential 
operator (or any other Green’s function if the overall region of interest is smaller than the entire space, and 
some specific far-field boundary conditions are involved.) The resulting construction, i.e., the corresponding 
collection of convolution-type boundary integrals, is often referred to as Calderon’s (surface) potential [1], in 
the particular formulation provided by Seeley [2], This potential is a combination of a conventional single- 
layer and double-layer potentials. The (vector-)function on the surface formed by the magnitudes of the 
corresponding singular layers, is called the density of the potential. 

In the classical potential theory for second-order differential operators, the single- and double-layer 
potentials (driven by scalar densities) are typically introduced and studied independently, for the sake of 
solving specific boundary-value problems, namely those of the Neumann and Dirichlet type, respectively. In 
the foregoing general framework, the “pure" single- and double-layer potentials are obtained as particular 
cases when the interface discontinuity has a special structure continuous function and discontinuous 
derivative 1 for the single layer (monopoles), and the other way around for the double layer (dipoles). This 
allows one to represent a given solution to the homogeneous equation on the specified domain as either only 
a single-layer potential, or only a double-layer potential, or, basically, any combination of the two. We note 
that different, representations of a given solution may have certain advantages with respect to one another 
for the applications related, e.g., to active control of time-harmonic wave fields, as explained in Section 4. 
In practice, to obtain the density of the single-layer potential one will need to solve a complementary 
problem of the Dirichlet type, and to obtain the density of the double-layer potential one will need to solve 
a complementary problem of the Neumann type, see Section 2. 

To implement similar considerations on the discrete level, we first need to introduce counterparts to the 
differential operator and its inverse, i.e., convolution with the fundamental solution (or a Green’s function). 
This is done as follows. Given a large region, we first mesh it and choose an appropriate finite-difference 
operator, and then select a far-field boundary condition at the outer boundary of the region so that the 
resulting finite-difference problem (a linear algebraic system) be uniquely solvable for any right-hand side. 
Next, analogously to the continuous case we define the discrete potential as a pair of independent grid 
functions, so that each one solves the linear homogeneous finite-difference equation with the operator that 
we have chosen on one of the two predetermined subdomains of the overall region of interest. We emphasize 

that, even when the subdomains have irregular shape, we can still use regular grids, which is convenient (see 
Section 3). 

Applying the finite-difference operator to the aforementioned pair of solutions throughout the overall 
region, we obviously generate a right-hand side on the grid. Specific values of this right-hand side, as 
well as specific grid subsets, on which it will differ from zero, will be determined by the stencil of the 
discrete operator and the behavior of each solution from the pair near the boundary that separates the 
subdomains. In Section 3, we will see that the latter behavior can be conveniently categorized in terms of 
“surface discontinuities” of the grid functions. As concerns the resulting right-hand side, for finite-difference 
operators of the second order it will always be defined on a two-layer “curvilinear” subset of grid nodes that 
follows the geometry of the continuous interface. This right-hand side is called the density of the discrete 

Hie direction of differentiation for the operators that we consider should be normal to the interface, see Section 2. 



potential. The potential itself is reconstructed by solving the foregoing finite-difference problem driven by the 
density. If, in so doing, we consider some special classes of surface discontinuities on the grid, see Section 3, 
then we can recover the discrete single-layer potential, for which the grid density is concentrated only on 
oik 1 layer of nodes, and the discrete double-layer potential, for which the density can be represented as a 
collection grid dipoles aligned with the coordinate directions and defined on a two-layer fringe of nodes. 
The key result in the discrete framework is that similarly to the continuous case, one can represent a given 
solution of the homogeneous finite-difference equation as a variety of discrete surface potentials, including a 
pure single-layer potential, a pure double-layer potential, as well as combinations of the two. 

Let us emphasize that our definitions of discrete surface potentials are not related in any respect to the 
notion of approximation of the continuous single- and double-layer potentials on the grid. Strictly speaking, 
we introduce the potentials of finite-difference operators rather than approximations to the potentials of 
differential operators. As such, our approach is fundamentally different from the previous techniques that 
have used mesh analogues of surface potentials, e.g., the method of capacitance matrices [3], in which specific 
boundary conditions (of either Dirichlet or Neumann type) are built into the construction of the grid potential 
from the very beginning, using interpolations between the regular grid and irregular continuous boundary. 
Nonetheless, the discrete potentials that we define here do approximate their continuous counterparts as well. 
To this end, we need to mention that the roots of our approach can be found in the work by Ryaben’kii [4], 
that has later developed into the difference potentials method (DPM), see [5,6]. It turns out that the discrete 
surface potentials the way we introduce 1 them in this paper can be considered a particular version of the 1 
so-called generalized difference potentials of [5,6]. As such, one can establish the approximation properties 
of the discrete surface potentials independently, using general apparatus of the DPM, sec [o <]. 

The motivation for defining surface potentials of finite-difference operators and analyzing their properties 
comes from the problems of active control of time-harmonic wave fields, in particular, common environmental 
noise, see, e.g., general reference [8], as well as our work [9,10]. In this framework, densities of the discrete 
surface potentials are conveniently interpreted as near-surface acoustic control sources that, are employed to 
cancel the adverse component of the overall acoustic field on a predetermined domain. In this perspective, 
depending on a particular setting, either monopole or dipole sources may be preferential from the standpoint 
of engineering implementation. Moreover, some of these source distributions may appear optimal with 
respect to particular criteria, e.g., the overall acoustic power or integral volume velocity, see [10]. 

The rest of the paper is structured as follows. A summary of results from the potential theory that 
are relevant to our current study is provided in Section 2. Discrete counterparts to the continuous surface 
potentials are introduced and analyzed in Section 3. Finally, in Section 4 w r e outline a mathematical frame- 
work for active noise control, and discuss some results in both continuous and discrete formulation, including 
the connection between the two formulations, which relates to the approximation properties of the discrete 
potentials. 

2. Surface Potentials of Differential Operators. In this section, we adopt the framework of dis- 
tributions as it will be convenient for introducing finite-difference analogues in Section 3. Further detail on 
classical potentials can be found in, e.g., [11,12]. 

Any second-order elliptic differential operator of the type 


Lu= 

>,j= 1 


d 2 u 

<llj dxjdx j 




( 2 . 1 ) 
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on 1" can be reduced to the canonical form 


Lu = Au + cu, (2.2) 

where A is the /(-dimensional Laplacian, by a non-degenerate linear transformation of the independent 
variables xj, j - 1, . . . ,n. Provided that the coefficients aij in (2.1) are constant, this transformation will 
not depend on the spatial location, 2 i.e., will be one and the same for all x € 1" . In this case, which is the 
only case that, we are going to study in this paper, we can effectively disregard the original form (2.1) and 
consider the operator L given by (2.2) from the very beginning, assuming, if necessary, that the corresponding 
transformation has already been applied. The coefficient c. in (2.2) is also assumed constant; the case c = —k~ 
(where k > 0) corresponds to the Yukawa operator, the case c = 0 to the Laplace operator, and the case 
c = k 2 to the Helmholtz operator. The dimension n of the space can be arbitrary, although from the 
standpoint of applications it is interesting to study either n = 2 or n = 3. The finite-difference analysis in 
Section 3 addresses primarily the case n — 2, for the reason of simplicity. 

Let us introduce the domain fl C R n and its complement 11, = D\fi to a larger domain ED C 1” , which 
nicu, in particular, coincide with the entire space R". il and 11, are the subdomains, and D is the overall 
large region, that we have referred to in Section 1. We also introduce the boundary between the subdomains 
and denote it. T = Oil To avoid possible uncertainties, we will assume that the domain ii is bounded, 
whereas its complement SI, . along with D, may be either bounded or unbounded. 

Next., we will need to introduce a special class of functions U defined on D, which will contain the 
solutions of the inhomogeneous differential equation Lit - f. The functions a £ U are assumed locally 
integrable on 0: i c L., 0 ‘(D) (as such, every u £ U is a regular distribution), although later we will 
restiic.t ourselves to an even more narrow set of functions as solutions, namely, piece-wise smooth, so that 
the operator L of (2.2) can be applied to u in the classical sense everywhere on D except, maybe, at the 
interface T. The functions u £ V shall also satisfy some linear homogeneous boundary condition on the 
outer boundary c® that would guarantee uniqueness of the solution, i.e., that the only solution of Lu = 0 
on ID), u e U, is trivial. We will sometimes refer to this boundary condition that essentially defines the class 
C as the jar-field boundary condition. The choice of the far-held boundary condition is typically not. unique 
for a given setting. Its specific form will in any event depend on the configuration of D and the type of the 
operator L of (2.2). In this connection we emphasize that whereas for any particular application we may 
need to consider the far-held boundary condition of an appropriate type determined by this application, for 
the general analysis of the current section we only need to ensure the uniqueness, and as such, any far-held 
boundary condition that provides it will be appropriate. For example, if D = E n , then for the Yukawa 
equation we only have to require that the solution vanish at infinity, for the Laplace equation it, either has 
to be bounded or, again, vanish, depending on whether n = 2 or n > 2. and for the Helmholtz equation 
is has to satisfy the so-called Sommerfold radiation boundary condition, see [11], (For that, the function 
u £ l has to possess additional regularity near 3D, namely, it. has to be at least (7* -smooth outside some 
large sphere, see [12].) If 0 is bounded, then, for example, the zero Dirichlet boundary condition at 3D will 
be appropriate, provided that in the Helmholtz case the domain D is not resonant, i.e., that —k~ is not an 
eigenvalue of the Laplacian on D. 

Let us first consider the unbounded case D = 1" . The fundamental solution of the operator L (see, 

-Ii there were first-order terms in (2.1) as well, they could have been eliminated by a slightly more elaborate transformation, 
see [11], so that the canonical form will still remain as given by (2.2). 
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e.g., [12]) is a solution £ = £(x) to the inhomogeneous equation 


L£ = S(x) 


(2.3) 


that is defined on the entire space R n and belongs to the corresponding class U: £ € U. For n = 3, the 
fundamental solutions £ G L r of the Yukawa, Laplace, and Helmholtz operators are given by the expressions 


£(x) = - 


-k\x I 


47t| rc | 


£(*) = 


47r|x| 


and £(x) = - 


o~ik\x\ 

4tt|x| 


(2.4) 


respectively. For n = 2, the corresponding expressions for c / 0 involve Bessel functions, and the fundamental 
solution of the Laplace operator is logarithmic, see [11,12]. 

For any distribution / such that the convolution £ * f exists in U , this convolution gives a unique solution 
u g U to the inhomogeneous equation Lu — f. Indeed, L(£ */) = L£ */ = £*/ = /; and uniqueness follows 
from the definition of class U. This, in particular, means that a given function u G U can be represented as 


u = £ * Lu , 


(2.5) 


provided, again, that the convolution exists in U. Indeed, denoting / — Lu and assuming that 3£ * / € T, 
we conclude, by the previous consideration, that this convolution, i.e., the right-hand side of (2.5), yields 
a unique solution of the equation Lu = / in the class U ; as such, it has to coincide with u. We also note 
that representation (2.5) is unique, i.e., if 3 g such that U 3 u == £ * g, then it has to be g — Lu. Indeed, 
subtracting the equality u = £ * g from (2.5) we obtain £ * (Lu - g) — 0, and then, applying the operator L 
we have L{£ * (Lu — g)) — L£ * (Lu — g) — S * (Lu — g) — Lu — g = 0. 

We will now focus on a particular subset of functions u = u(x) G U , x G M n , namely: 


w(x) 


t;(x) for x G Q 
w(x) for x G Hi 


( 2 . 6 ) 


where Lv - 0 on H, and Lw = 0 on In other words, u(x) is composed of two independent branches 
v(x) and ?e(x), such that u(x) solves the homogeneous equation on the domain ft, and w ( x ) solves the 
homogeneous equation on its complement ftj and satisfies the appropriate far- field boundary condition. On 
the interface T the function u(x) is, generally speaking, discontinuous. 

It is clear that the application of the operator L of (2.2) to the function u of (2.6) yields a singular 
distribution with the support on the interface T only: supp Lu C F. Specifically, let (w, (j>) denote the linear 
functional associated with a given distribution u, here </> is a test function. Then, differentiating u of (2.6) 
once in the sense of distributions, we obtain 


(V«,0) = -(u,V0) = - [ u(x)V<f>(x)dx = - v(x)V <j)(x)dx - / w(x)V<f>(x)dx 

J]R n «/^ 2 JQ\ 

= - f V (v(x)<f>{x)) dx + [ <t>{x)Vv{x)dx - / X? (w(x)4>{x)) dx + / 0(x)Vu>(*)d* 

Jo Jn Joi J^! 

= / 0(x) (Vu(x)} dx - / ?;(x)4>(x)nd<r + / uj(® )<£(«) nd<7 

Jj Jr 

= [ 0(x) {Vu(x)}dx + f [u]r(x)(/>(x)nda = ({Vu(x)} 4* [u] r n<$(r), </>) , 

Jr« Jr 


where dx and da correspond to the volume and surface integration, respectively, n is a unit normal to T 
pointing toward fb, {Vu(x)} denotes the regular part of the gradient calculated in the classical sense where 
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is exists, i.e., on Q U fi u and [u]y stands for the discontinuity of the function u{x) across the interface: 
Mr(®) — lirn^^j. te(x ) — lim a .^_ )>x v ( x ), where x 6 T, x f £ fij, and x tf £ fi. In other words, we have 

Vu = { Vu(x) ) + [u] r n<5(r). (2.7) 

For every coordinate direction Xj, j = 1, . . . , «, the second term on the right-hand side of (2.7) is a single 
layer on the surface T with the magnitude [u]r(x)?ij (x), x £ F. Let us now apply the divergence operator to 
both sides of equality (2./), and when differentiating the function {V'u(x)}, which is possibly discontinuous 
across T, implement the same rules and notations as before. This yields: 


Au = {Au} + <[Vtt] r , n)6(T) + (V, [ u ] r nS( r)>, (2.8) 

where (*, *) denotes a standard dot product of n-dimensional vectors. The second and third terms on the 
light-hand side of equality (2.8) can be rewritten using conventional definitions of differential operators in 
the space of distributions: 


({[yu]r. n)rf(T), 0) — J ([Vu]r, n)(bda = J ((Vw,n) — (Vr. n))<f)da 


dw Ov 
y \dn On, 


du 

dn 


m,<i> 


where normal derivatives and on T shall be understood as uniform limits from the side of fl } and H, 
respectively, provided F is sufficiently smooth, and 


da 


((V,[u}yn6(Dl<t>) = “il (m r ”^(r)’^) = - iz J r ^ rnj ^j 

= - Jw\ r(V0, n)d(T = - j[u\Y^da = (^([«]i ^(F)), <£ ) . 


Altogether, from equality (2.8) we obtain 


d 

HD + ^(Mr<5(r)). (2.9) 

The second term on the right-hand side of (2.9) is a single layer on the surface T with the magnitude 
[s^Jr x e r, and the third term on the right-hand side of (2.9) is a double layer on the surface T with 
the magnitude [u] r (x), x € T. Finally, for the operator L of (2.2) we have 


Au = { A?i} + 


du 
d n 


Lu = 


du 

On r 


S(T) + 


d 

dn 


U«H(r)), 


( 2 . 10 ) 


because the zeroth-order term cu contributes only to the regular part, and {Lu} = 0. 

The right-hand side of expression (2.10) is compactly supported; consequently, its convolution with the 
fundamental solution £ exists in the sense of distributions. Moreover, it can be shown that this convolution 
does belong t.o the class U (see [11, 12]), i.e., satisfies the far-field boundary condition at infinity that is 
appropriate for every particular type of operator (2.2). As such, we can use formula (2.5) and obtain 


u = 



du 

Ihi 


(y) 


M r (y)|^(z - y)) 


da TJ 


( 2 . 11 ) 


The integral on the right-hand side of (2.11) reconstructs the original function u(x), i.e., it coincides with 
v(x) for x £ {} and with ?e(x) for x £ This integral is called the surface potential , and it has two 
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components — a single-layer potential and a double-layer potential. Hereafter, we will, for simplicity, refer 
to the vector function ([f^] r , Mr) 7 011 the interface T as to the density of the potential. We reiterate, 
however, that the actual density of the potential is a. singular distribution [^] S(T) + 777^ ( [t/]<5(F)) with the 
support on F, which, by its nature, is a right-hand side to the differential equation that we study. 

We see that the density of the potential is fully determined by the discontinuity of the function u of 
(2.6) itself and its normal derivative. This is a key consideration in our analysis, namely that the behavior 
of a piece-wise smooth solution to the homogeneous equation is completely controlled by the location and 
type of discontinuities that it has. A similar consideration will be employed in Section 3, when discussing 
the discrete formulation. We also notice that when w(x) = 0, x G Hi, then formula (2.11) reduces to the 
standard Green’s formula that recovers the solution u of the homogeneous equation Lu = 0 on the domain 
ft through its values on the boundary T. Reciprocally, if v{x) = 6, x e ft. then (2.11) becomes the Green’s 
formula for the complement domain fti. 

;From the definition (2.11) one can see that changing the function w(x) on the domain fti does not 
change the values of the potential on the domain ft. As such, the following natural question arises: Given a 
fixed function e(x), Lv = 0. x G ft, describe all its possible representations on ft in the form of the potential 

v = Jr (s(x - y)Zo(y) - da, J' ( 212 ) 

where £0 and £1 are scalar densities defined on the boundary T. The answer to this question is given by 
Proposition 2.1. Surface integral (2.12) coincides with the given function v(x), Lv — 0, on the domain 
ft, if and only if £ 0 = §^7 - and £1 = w - v on T, where w{x) is some solution to the equation Lw = 0 

that is defined on fti and satisfies the corresponding far- field boundary condition. 

In other words, by considering the entire variety of appropriate auxiliary functions w(x), x G fti , we 
obtain all possible representations of the given u(x), x e ft, in the form of the potential (2.12). 

Proof. One implication has already been established, namely, we have constructed the potential (2.11) 
so that it reconstructs the given v for any w, Lw = 0 on fti . Conversely, assume that v is represented by 
formula (2.12) on ft. We need to show that there is a function w that solves equation Lw = 0 on fti, satisfies 
the far-field boundary condition, and such that £1 and £0 will be discontinuities between w and v on the 
surface F ; more precisely, £1 will be discontinuity of the function itself, and £0 will be discontinuity of the 
normal derivative. Consider the standard Green’s formula for u(x), x € ft, 

V = | r (-£(* - - »>) *». 

and subtract it from the general representation (2.12). The result will be identically equal to zero on ft. On 
the complementary domain fti, the resulting function 

w = J U(x - y) ^o(») + f^(3/)) - (6(y) +u(»))|^(a: - y)j da v (2.13) 

satisfies Lw = 0 because the source terms are concentrated only on T. The far-field boundary condition can 
again be verified directly for every specific case (\ukawa, Laplace, Helmholtz) by analyzing properties of 
the integrals, [11,12]. Consequently, the foregoing general argument on representing a given function in the 
form (2.5), or more precisely (2.11), can be applied to the function that is identically equal to zero on ft and 
equal to w of (2.13) on fti, which allows us to conclude that £o + f” = an(1 £i + v = w 011 r - D 

As Proposition 2.1 describes the entire variety of representations for v in the form of a boundary potential, 
one can easily identify some important special cases. For example, we can represent v(x), x € ft, as a single- 
layer potential only. To do that, we need to find w(x), x £ fti, such that the overall function u(x) of (2.6) 
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will have the discontinuity only in its normal derivative, and not in the function itself. This w(x) will be a 
solution of the following external Dirichlet problem: 

Lw = 0, x 6 ft,, 

I , r . (2.14) 

w 'lr = w li” w e b, 

where the inclusion w € U simply means that w has to satisfy the appropriate far-field boundary condition. 
Problem (2.14) is always uniquely solvable on ft, = !"\ft. Similarly, v(x), x e ft, can be represented as a 
double-layer potential only. In this case, the function w(x), a: € ft,, has to be such that u(x) given by (2.6) 
will have discontinuity only in the function itself, and not in the derivative. This function w(x) shall solve 
the following external Neumann problem: 


Lw — 0, x £ | , 

dw du f (2.15) 

, weu. 

on r on lr 

We note* that the additional necessary solvability condition of zero total flux through the interface that needs 
to be imposed for the two-dimensional Laplace equation in problem (2.15) is satisfied because the boundary 
data themselves come from the harmonic function on ft. Of course, besides the two “polar” cases of a 
pure single-layer potential and pure double-layer potential, the function v(x) on ft can be represented as a 
variety of combinations of the two potentials, see Proposition 2.1. As described in Section 4 below, different 
representations of a given field v in the form of a surface potential may exhibit different properties in the 
framework of active control of sound. 

Let us also note that the argument used when proving Proposition 2.1 provides for another useful way 
of looking at the properties of surface potentials. Consider an arbitrary pair of functions £ 0 and £ on T and 
introduce £ = £o<i(T) + ^(£]<5(T)). Then the potential with the density £ is given by u = £ * or by the 
integral on the right-hand side of (2.12). Obviously, it satisfies the equation Lu = £, and one can also show 
that a e l . Then, £ 0 = [f^] r and 0 = [«]r If we now denote v(x) = u(x)| i€S; and w(x) = u(x)| i€n as 
in (2.6) then we conclude that for a given £ on T the potential £ * £ renders die decomposition of | into a 
part ( 1 1 " ' ll ) that can be extended to ft so that the extension u(x) satisfies Lv = 0 on ft, and the 

lemaining part ( ^ |,., w|r) that can be extended to ft, so that the extension w(x) satisfies Lw = 0 on ft. 
In the framework of the Helmholtz equation the meaning of the field v is the incoming wave with respect to 
the domain ft, because Lv = 0 on ft, and as such v can be attributed to some (unknown) sources outside ft, 
i.e., on ft, . Reciprocally, w shall be understood as the outgoing wave with respect to ft. We therefore can 
sa> that for a given £| r that consists of the traces of both incoming and outgoing waves, the potential £ * £ 
automatically reconstructs only the incoming part of the field on the domain ft and only the outgoing part 
of the field on the domain ft, . This capability is important for active noise control, as explained in Section 4. 

To conclude this section we should only notice that the case of a bounded domain D can be analyzed 
similarly. Instead of the fundamental solution £ we will need to use the Green’s function Q that satisfies 
a specific far-field boundary condition on 5D (for example, the zero Dirichlet boundary condition). The 
Green s function is a function of two (n- dimensional) variables x and y, Q = Q(x, y)\ it. is defined as 

&( x , y) — ~ - y) + X(x, y), VyeD: G{x,y)£U, (2.16) 

where Vj, e D : L x \(x,y) = 0, x € ©. In other words, the Green’s function is obtained by adding to 
— c-(x — y ) another function Xtx. y ). which is a solution to the homogeneous equation with respect to the 
variable x on D, and such that the sum satisfies the boundary condition on <91D> with respect to x for any y. 
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When one uses the Green’s function, the expressions for surface potentials in the form of boundary integrals 
remain the same as before with only —£(x — y replaced by Q(x, y): and the result of Proposition 2.1 also 
holds. Of course, in this case one needs to remember that every time a function w on is considered that 
satisfies the far-field boundary condition, this boundary condition is a boundary condition at <9B rather than 
at infinity. As far as uniqueness of the solution for problems (2.14) and (2.15) in the new space U, it follows 
from the maximum principle for the Yukawa and Laplace equations, and for the Helmholtz equation we need 
to assume that the complementary domain 11] is not resonant. 

The use of the Green’s function Q of (2.16) as an alternative to the fundamental solution £ allows 
one to take into account specific, boundary conditions on <9D, in case this is important. As concerns the 
representation of v(x) on G in the form of a surface potential, replacing £ by Q does not change the results 
in any respect. Consequently, one can use a specific Green’s function instead of the fundamental solution 
not only when it is necessary, but also when it is convenient. In particular, in the following Section 3, when 
discussing the discrete formulation, we adopt the framework of a finite domain ID) with the zero Dirichlet 
boundary condition on its boundary 9B. 


3. Surface Potentials of Finite-Difference Operators. Let us introduce a discrete grid N on the 
domain D. For simplicity we can assume that D is a large rectangle aligned with the Cartesian coordinate 
directions, and N is a regular square-cell grid with the size h. We will denote by u^'K etc., the 

discrete functions on the grid N. Having constructed the grid and defined the grid functions, we discretize 
the differential operator X of (2.2) and denote its discrete counterpart by L [h) . In this paper, we use 
finite differences to construct discretizations, and consider only second-order finite-difference operators. The 
simplest of the latter involves standard central differences for the Laplacian, and in the two-dimensional case 
is built on the symmetric five-node stencil: 


tf h) u w 


i {h) 

l i±Uj 


+ U 


(h) 

i-lj 


4- u 


{h) 
ij ± 1 


+ U 


(h) 

u-i 


_ 4u ih) 


h 2 


+ n b‘j- 


(3.1) 


We emphasize that in general we can consider more complex grids and domains, as well as more sophisticated 
discretizations, but for the purpose of demonstrating the concept, hereafter we rather adhere to the most 
straightforward constructions, namely, a uniform Cartesian grid on a rectangle ID) and the central-difference 
discretization (3.1). Moreover, the apparatus that we develop in this paper, as well as the general apparatus 
of difference potentials [4 7], provide a viable means for dealing with complex geometries on regular grids. 
As such, the discrete formulations that we analyze in this section should not be regarded as model examples 
only, but rather as an outline of the direct finite-difference approach to analyzing the applied problems (see 
Section 4 for more detail). 

Similarly to the continuous case, the operator L {h) needs to be supplemented by the far-field boundary 
condition on cW. It is natural to assume that the outermost collection of nodes of the grid N belongs to the 
rectangular perimeter OID). Then, imposing, for example, a zero Dirichlet boundary condition is eas\. If we 
choose a generic notation n for the nodes of the grid N (n, in fact, is a pair of indexes (i, j)), then we simph 
say that 


u ( n h) = 0 for all n€{Nn<9D}. 


(3.2) 


In general, we define the space of grid functions U^ h) as all those and only those functions on N that satisfy 
the discrete boundary condition on 3D that we have selected. 

Besides the grid N and the space of functions U {h) defined on it, we will also need to introduce the grid, 
on which the residuals of the operator of (3.1) are defined. Subsequently, the discrete right-hand sides 
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of the inhomogeneous finite-difference equation will be considered on this grid. We will call this new grid M. 
and the generic notation for its nodes will be m (m is also a pair of indexes). The structure of the stencil 
of L {h - . see (3.1). implies that ni € M if and only if 3n € N: n 9© and m = n. In other words, for our 
simple formulation the grid M can be obtained from N by throwing away the outermost nodes of N that 
belong to 9®. We assume that for any discrete right-hand side f [h) = f { ,!‘ ’ , m € M. the inhomogeneous 
finite-difference equation £ ,A) ti (h> = f {h) is uniquely solvable in the class U {h) . If the class U (h) is defined by 
the zeio Dirichlet boundary condition (3.2) (this case is our primary example), then the unique solvability 
for the \ukawa and Laplace equations follows from the maximum principle. As for the Helmholtz equation, 
we need to additionally assume that the domain is not resonant, as we did in the continuous case. The 
aforementioned unique solvability allows us to introduce the inverse operator to L (h) . We will denote it G (h) 
so that for a given discrete right-hand side / (/l) = m e M, the solution U {h) 3 u (h) = u ( ,!‘ > , n e N, 
of the finite-difference boundary-value problem on N is given by u (/l) = G (h) f [h) . The operator G [h) is a 
finite-difference analogue of the convolution with the Green’s function Q described in Section 2. 

Let us now introduce a special notation for the stencil of the difference operator of (3.1). For 
m = (i,j, ) the stencil N,„ centered at m e M consists of the five nodes of the grid N: N,„ = {(i,j),(i + 
1 - 1)}- Next, in accordance with the partition of domain © into 11 and its 

complement Jli, we define the corresponding subsets of the previously constructed grids: 


M + =Mnn, mt = m\m + = m n n, , 
N + = U Nm, N- = U Nm, 

me wr 

y = Ttr\N~ J = N“ n n, 7" = N + n Q, . 


(3.3) 


We emphasize, that the grid M that pertains to the right-hand side of the finite-difference equation is 
partitioned into M + and M“ directly, i.e., following the partition of D. In contradistinction to that, the grid 
N is not partitioned directly, we rather consider the collection of all nodes of N sw f ept by the stencil N w when 
its center belongs to M+ , and call this sub-grid N+, see (3.3). Obviously, some of the nodes of N+ obtained 
this wa\ are outside Q, i.e., in Qi, and these nodes are called y — . The sets and 7+ arc 1 defined in a 
similar way starting from . The key idea is that whereas the grids M+ and M” do not overlap, the grids 
N and N do overlap, and their overlap is denoted 7 i obviously, 7 — 7^ U 7~ . The subset of grid nodes 7 
is c alled the grid boundary, it is a two- layer fringe of nodes'* that is located near the continuous boundary T 
and in some sense straddles it. The geometry of the domains and the grid boundary is schematically shown 
in Figure 3.1. 

Next, we consider a grid function u w = Un l) € U {h) 


u 


ih) 

n 



for n e N + 
for n € N“ 


(3.4a) 


where = 0 on M + and L^ h \v^ h) — 0 on M“ . Similarly to the c*ontinuous case, see (2.6), the function 

u {fl) of (3.4a) is composed of two independent branches on N + and N~ , so that each of them solves the 
homogeneous finite-difference equation on the corresponding sub-grid M + or M“ , respectively. O11 the grid 
boundary 7 the function u {h) is defined as double- valued. The notion u {h) e U {h) only means that w {h) 
satisfies the far-field boundary condition that we have chosen, e.g., (3.2). 

3 For more elaborate stencils it may be a “thicker” multi-layer fringe. 
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Fic;. 3.1. Schematic geometry of the domains and the grid boundary: hollow bullets solid bullets - 7 • 


Along with considering the function of (3.4a), we will also need to eliminate the ambiguity in its 
definition on 7 and introduce = u\i * € U^: 



-(h) 

Vn 

for 

n £ N+\ 7 “ 

-(h) 

Wn 

for 

n £ N - \ 7 + 


(3.4b) 


where and w {h) are again solutions of the homogeneous finite- difference equation, but considered oil 
smaller grid subsets; namely, X^fd^ — 0 on M + \ 7 + and = 0 on M \7 . The function u * * of 

(3.4b) is uniquely defined everywhere. We emphasize though, that formula (3.4b) defines a class of functions 
which is wider in some sense than that defined by (3.4a). Indeed, every function u ^ ^ of (3.4a) can be 
reduced to the corresponding u {h) of (3.4b) simply by means of truncation; however, given a function u (h) of 
(3.4b) we will not necessarily always be able to extend its branches v {h) and w {h) from N + \ 7 ~ and N"\7 + , 
respectively, to N + and N“ , respectively, so that the extensions satisfy the homogeneous equation in the 
same sense as the branches of do. This is easy to see already from the following simple example. Assume 
that we have a solution to the homogeneous difference Laplace equation that we want to extend from N \q + 
to N“. Even if this particular extension does exist, we cannot, generally speaking, continue extending it 
further and further “inwards. ’ Otherwise, we will eventually end up with a solution to the homogeneous 
equation on the entire grid N that satisfies boundary condition (3.2), which will either violate the uniqueness 
(or maximum principle) or should be trivial. 

The single- valued function u {h) of (3.4b) will be needed mostly for “technical” purposes, namely, to first 
apply the finite-difference operator L {h \ and then use the inverse operator G (/l) , which altogether will allow 
us to define grid analogues of continuous surface potentials. As for the function of (3.4a), its double- 
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valued property on 7 can be conveniently termed as “discontinuity” of the grid function on the grid boundary 
7. In this perspective, the function of (3.4a) will allow us to analyze various kinds of such discontinuities 


and accordingly identify different types of discrete surface potentials, in particular, the single-layer potential 
and the double-layer potential. 

Applying the operator £<*> of (3.1) to the grid function u {h) of (3.4b), we produce the right-hand side on 


the grid. By construction, this right-hand side may only be concentrated on the grid boundary 7. Indeed, 
consider, for example, the interior sub-grid M + . For every m G M+ such that N w D 7“ =0, we obviously 


. and v {h) is a solution of the 


have L {h) u {h) \ = 0, because for all such m’ s clearly L {h) vS h) = L {h) v [h) 

homogeneous equation on M + \7+. Consequently, the only nodes m € M + , for which L {h Ht {h) J may differ 
from zero, are those that satisfy N m Pi 7“ ^ 0. By construction of the grid sets (see Figure 3.1). these are 
the nodes 7 - Reciprocally, the only nodes 111 G M , for which may differ from zero, are those 

that satisfy N m H 7+ ^ 0, and these are the nodes 7“. Altogether, we can write 


L {fl) u {k) — f\ h \ where Vm G 


fit *> 


= 0. 


(3.5) 


Next, applying the inverse operator G {h) to the grid right-hand side f\ tl) of (3.5), wo obviously restore 
the original function u ( h ] on N (due to the aforementioned unique solvability of the discrete' problem): 


fiW _ Q{h) jr(h) 


( 3 . 6 ) 


We will call the representation of u {h) in the form (3.6), where the source term f (h) is concentrated only on 
the grid boundary 7? the discrete surface potential with the density /i/ 2 2 . Similarly to the continuous case, 
see formula (2.11) and subsequent discussion, the density of the potential by its nature is the right-hand side 
to the equation. 

Now in the discrete framework, we are ready to formulate a question similar to that previously addressed 
by Proposition 2.1 in the continuous case. Namely, for a given function on N" 1 "\7~, £i h )^(h) q Qn 
M + \7 + , describe 1 all its possible representations in the form (3.6), where again, the source terms can only 


differ from zero on the gird boundary 7. The answer to this question is given by 

PROPOSITION 3.1. Expression (3.6), where supp f\ h] = 7, coincides with the given function v 


v ^ Jh) f 

L w V h) = 0, on the grid [f\f, if and only if f\ h) = L w u {h) , where «<*> is composed of v ih) and 
according to definition (J.Jb), and ?/’ ' ^ 1 is a grid function on N _ \7~*~ that satisfies the equation — q 

on M~\7~ and the selected far -field boundary condition, e.g., (3.2). 

Proof. One implication has already been shown, namely that for any given ^ and arbitrary 
the potential (3.6) wdth the density f\ h ^ obtained through (3.5) will reconstruct the original on N + \7 - . 


Conversely, assume that there is / 


•(h) 


on 


such that fy h) = 0 on M\y, and G {h) f 


>( h ) 


N+\v 


= v^ h - . Consider 


u 


(h) 

n 


v\i ] for n£N + \7 
0 for n G \7+ 


(3.7) 


and define 
have I 


u [h) 


= r h} u [h) ; clearly, g\ h) = 0 on M\y. Let w ( '*> = G ih) (/' ft) - Because of (3.7), we 

— 0. Moreover, = 0 on M \y since both f\ 11 1 and g ^‘ 1 may only differ from zero 


on 7. Finally, w 1 ^ 1 satisfies the far-field boundary condition by construction. If we now assemble u ' ^ ■ of 

and w ill] according to (3.4b), then it is easy to see that u [h) - u (h) + w (h) on the entire N. Consequently. 
L lh) u< h) = L (h) fl (h) + L (h) G (fc) _ g W\ = g W + jU,) _ ; th) _ Q 
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iFroin Proposition 3.1 we can conclude that by varying the auxiliary function w w on N \7 + we obtain 
all possible representations of the given on N + \7“ in the form (3.6), where supp f\ h) — 7. 

Let us now return to the more narrow class of functions ) defined bv (3.4a): these functions are double- 

valued on 7. First of all, using Proposition 3.1 we can easily see that for a given function defined on 
N+ such that L {h) v {h) = 0 on M+ , its truncated portion t? (/l) | N+v> _ can be represented as a discrete surface 
potential (3.6) with the density }\ h) concentrated on 7, and all such representations are parameterized by 
the function w {h) defined on N“\7 + . Formally, we can say that the entire v {h) on N + is represented as a 
discrete surface potential; for that, we only need to supplement f\ h) with the values _ that 

N + \7“ 

are known. Next, we recall that w ^ of (3.4a) is, in fact, defined on N _ . This will allow us to analyze the 
structure of the right-hand side f\ h) of (3.5), i.e., the density of the potential. 



1 1 1 r 


Fig. 3.2. Schematic for a fragment of the grid boundary: hollow bullets — 7 + , solid bullets 7 . 


Consider a fragment of the grid boundary schematically shown on Figure 3.2 and take an arbitrary node, 
say, node b . If we now denote the truncated functions: and = u f(,l) | N _^ +T and have 

u {h] defined according to (3.4b), then we obviously obtain 


ih) (h) (h) (h) 

/(Ml =L WfM _ LW V U>)\ + Wa + w, ‘ , , 


wiV-ifr vP-vM [«,<*> -v‘*> vP-v™ 

h 2 h 2 K 2 h 2 
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Similarly, for the neighboring node c we have 


fl 


(h) 


(h) 

Wc 


Jh) 


Jh) 


Jh) 


h 2 


h 2 


{h) 

nr c - V , 

T 2 


ih) 


W 


\ h) _ ,JA) 


h 2 


(3.9) 


W 


at 


;From formulae (3.8) and (3.9) we first conclude, that for every given node from 7 , the value of f, a , 
this node gets contributions from those and only those neighboring nodes that are on the other side of 
the continuous interface T. More precisely, if m € 7 + , then the value of f\ h) | will be affected by the 
values of v {h} and w {h) at those and only those nodes of y~, for which N„, n y~ ± 0, as well as bv «<*>| 
itself. Reciprocally, if m € y~ , then the value of will be affected by the values of w <*> and t>‘*> at 

those and only those nodes of 7 +, for which N,„ D 7 + / 0, as well as by itself. In other words, 

the “homogeneity is generated when the finite-difference stencil spans across two different solutions of the 
homogeneous equation. It is generated only on 7 because everywhere else on the grid the stencil applies to 
one branch of «<*>, either or «;<">, and does not. touch the other one. In this sense, it will be convenient 
to say that the discrete function fi ih) of (3.4b) is “ discontinuous ” on the gnd boundary 7 . And the double- 

rained property of its parent function on 7 (see (3.4a) will, in fact, help us identify different types of 
such discontinuities. 

To do that let us first consider the segment of the grid line between the nodes b and e, see Figure 3.2. 
Looking at formulae (3.8) and (3.9), one can say that this segment contributes both to the value of /i/' ; 

and to that of f\ *| . Indeed, the contribution of the segment [6, c] to /- 


.(/>) 


IS 


fi h) 


(M) = 


wl h) -vl h) 


.(h) 




h 2 


and the contribution of the segment [/;, r] to f- 


( h ) 


(3.10) 


A h) «M) = 


w< h> - vi h) 

h? 


(h) 

Wr 




h 2 


(3.11) 


To obtain the full value of /i/ 0 ^ one obviously needs to add up the contributions from all segments that 
originate in b and intersect T , i.e., = f\ h (M) + /<'° [([&,«])■ Similarly, the full value of f \ h) | is 

given by all segments that originate in c and intersect T, i.e., f\ h) = f\ h) ([6,c]) -i- fi h) ([d,c]). and the 
same is, of course, true for every node from 7 . This is just another way to 'say that for a given nod< 


ill "y. 


the value of f\ h) is affected by the neighboring nodes on the other side of the interface. Alternatively, one 
can say that all those and only those grid segments that have one endpoint in y+ and the other one in 7 ~ 
contribute toward the values of the right-hand side f\ h \ 

For a given grid segment of the foregoing type (i.e., with one endpoint in 7 + and the other one in 7 

setnnent ffr.rl. i t m 'A V in nnrtienlor 4 - ^ *- *(A)| /fL cf/l) I /r, i\ ^ 


say, segment [6,c], it may happen, in particular, that /• 

(3.11) show that, a necessary and sufficient condition for this is given by the equality 


([&,c]) — -f-, h) ([&,c]). Formulae (3.10) and 

b 


W h Hv)- v ™- l,W 

h b <- (t) ~ v 


b _ 


(h) 

IV C 


■ W 


h 2 


} (h) 

- 


(3.12) 


In other words, if we have a grid segment with the endpoints on different, sides of T, and if we want the 
contributions of this segment toward the values of /‘ fc> at both endpoints be equal in magnitude and have 


14 


opposite sign, then we have to require that the difference quotients Vj^(u) and V{^(w) of the original 
functions v {h) and w {h) along this segment be equal, V^(i;) = V^(u>) = V^, see formula (3.12). In this 
case, it is natural to associate with the segment [b,r] a discrete dipole with the moment (p, r : 


® c=ft '[ / ‘i fc) L([ 6 ,c ])" / ‘ fc) L([ 6,c ])] = 


n\ 


(h) 


— V 


(h) 


- V 


(h)\ 


(3.13) 


For the particular example that we have analyzed, the dipole (p )r of (3.13) is aligned with the Cartesian 
direction x , see Figure 3.2. Obviously, the construction of a dipole aligned with the direction y would have 
been the same. 

Assume now that for every segment with endpoints in 7 + and 7“ the corresponding difference quotients 
of v {h) and w (h) are equal. In other words, 


V (u 1,712) such that ?ii € 7 + , nj 6 7 fl N ni : 
.(h) Jh ) (h) ih) 


1 in' 1 — in' ’ 1 

IVTih ) / x _ l n } _ W112 ^ _ £ v (/*) / x 

Then, we can associate a pure discrete dipole with every such segment [711,7*2]: 

(ft) (/O 


(3.14) 


(In 1 n 2 — b 

• /n 

1 772 


_ rylh) 



— V m 77 2 V ( 

) — V U] 71 2 


([«1,» 2 ]) 

= -2 

71 1 



- V 


(3.15) 


of f~ h 1 from dipolo moments by tlie formulae 


V j 6 7 + : — E — 


1 n 1 2 /i 


V n*2 G 7 : /; 


(ft) 


E - 

^ 9/1 


(3.16) 


71 2 2/l 

7i 1 G7 + nN ri2 


1 712 ■ 


Thus, we can give the following 

Definition 3.1. If the functions v^ h) and u^ h) of (3.4a) satisfy boundary conditions (3.14), while the 
function of (3.4b) is obtained simply from the truncated components: 7 > ( h 1 = and — 

then the corresponding f[ h ^ of (3.5) is called the density of the discrete double-layer potential , 
and expression (3.6) defines the discrete double-layer potential itself. 

As we have seen, once the boundary conditions (3.14) hold we can equivalently re-define the density 
of the discrete double-layer potential through a collection of dipoles f/ ni „ 2 on the grid boundary 7. These 
dipoles (3.15) are associated with the grid segments that have one endpoint in 7 + and the other one in 7”. 
In contradistinction to the continuous case of Section 2, in which the orientation of the dipole layer was 
normal to the boundary T, here we work on a specific grid, and the discrete dipoles are aligned with the 
Cartesian coordinate directions. The point-wise values of the potential density can be reconstructed from 
the dipole moments using formulae (3.16). Let us also note that hereafter it is important to maintain that 
n\ € 7 + is the interior endpoint, and 112 € 7~ is the exterior endpoint of a given segment [7*1,712]; this way 
one can easily verify that all the difference quotients that we are considering, see formulae (3.14)-(3.16), will 
always be taken with the correct sign, no matter whether it is actually a forward difference or a backward 
difference in every particular instance. 
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Definition 3.1 is a natural finite-difference analogue to the notion of the* continuous double-layer potential 
introduced in Section 2. Indeed, a pure continuous double-layer potential is defined by the discontinuity [«] r 
in the function itself, and not in its derivative. Accordingly, boundary conditions (3.14) that Definition 3.1 
hinges upon essentially say that there is no derivative- type discontinuity in the grid function u {h) on the grid 
boundary 7 because the first difference quotients on all segments from 7 are equal from both sides. On the 
other hand, conditions (3.14) do allow, generally speaking, for a discontinuity of the function itself, and this 
is what may generate a non- trivial density of the potential in this case. 

LEMMA 3.2. Any v {tl) defined on N + such that L ih) v {h] — 0 on M+ , where L {h) is the operator of (3.1) 
with v < 0, i.e, either a Yukawa or Laplace operator ; can be represented in the form of a discrete double-layer 
potential. 

Proof What we essentially need to show is that for any v {h) that meets the requirements of the Lemma, 
there will always he a function w {h) on N - that solves the homogeneous equation L {h) w {h) = 0 on M“ , and 
satisfies the far- field boundary condition (3.2) and boundary conditions (3.14) on 7. The problem of finding 
the appropriate function w w on N - can be called a discrete exterior problem of Neumann’s type, because 
boundary conditions (3.14) are formulated for the first difference quotients. In this sense, this problem is 
analogous to the continuous problem (2.15) that needs to be solved in order to represent a continuous v(x) 
on il in the form of a double-layer potential. 

To show that the discrete problem for on is always uniquely solvable, we will consider the 
corresponding homogeneous problem 

L {h) w {h) = 0 on M“, 

(3 IT) 

v (»i, »•.>), »2 € 7 «1 € 7 + n N n , : wil‘J - «?<*' =0, 

and prove that its only solution is trivial. First of all, it is easy to see that the number of equations and 
the number of unknowns in problem (3.17) are the same; they are equal to the number of nodes of the 
grid M . Next, to show that the system is non-singular, we are going to employ a maximum principle 
argument. Assume that, there is a noil-trivial solution w (/,) to problem (3.17) on the grid and denote 
u- 1 = | N -\-, + - The truncated function ir ' ' can reach neither its positive maximum nor its negative 

minimum at any interior point of the grid, on which it is defined, i.e., at any node from M \7 — . We will 
prove the case of a positive maximum, the other case can be analyzed similarly. Take some rn € M“ Y)~ 
such that B n € N m \{m} : |^ . Obviously, if no such m exists, then the value of the function 

w (/l) is constant on the entire grid, and because of the far- field boundary condition (3.2) this constant is 
equal to zero. On the other hand, having found such a node m £ M - \7~ with positive maximum, we apply 
the operator L {h) to the function w w at this node and immediately see that the result can only be negative, 

| ^ <0, which contradicts one of our previous assumptions, see (3.17). Consequently, may only 
reach its positive maximum on the boundary 7“, and moreover, if w {h) £ const, and m € 7* is one of the 
boundary nodes at which the maximum is reached, then V n 6 N m n W - } is has to be I < tj&W i 

1 x * In \m' 

(Otherwise, it will be an interior node with the maximum value.) Calculating the value of L {h) w {h) j we 
obtain zero contributions from all those segments of the stencil N m that have the other endpoint (besides m 
itself) in 7 + because of the boundary condition on 7 in (3.17); we may obtain zero or negative contributions 
from all those segments (if any) that have the other endpoint in 7“ because is a maximum; and the 
contributions from all the segments with the other endpoint in MT\7“ will be negative. As such we again 
arrive at a contradiction: L {h hl^ h) | fji < 0, which proves the Lemma. □ 

To prove the result, similar to that of Lemma 3.2 for the Helmholtz equation, one actually needs to 


16 



assume that the corresponding exterior Neumann-type problem is solvable (i.e., there are no resonances of 
the complementary domain); then, the justification is straightforward. 

Let us now return to formulae (3.10), (3.11) and assume that condition (3.12) does not hold. This means 
that we will not be able to associate a pure discrete dipole with the grid segment [ft,c], because we will no 


longer have the property f\ h } ([b,c]) — -f\ h) ([6,c])- Instead, we can write 


A h) 


(M) = 

= ~A 


f(h) 

Jl 

(h) 


(M) + 


{ h ) 
W c 




(h) 


Jh) 


Aft) 


h 2 


fi 2 


(3.18) 


(M) + p| c ([M) 


The term p| ^ ( [fc, o] ) in (3.18) accounts for the discrepancy between the first difference quotients of the 
functions and along the grid segment [6, r]. This is analogous to the discontinuity of the first 
derivative (with respect to x , see Figure 3.2) in the continuous case. Recall, in the continuous framework 
(see Section 2) the discontinuity in the function [u ] r gave rise to the double layer, whereas the discontinuity 
in the derivative [^] gave rise to the single layer. Similarly in the discrete case, we can still associate a 
dipole qi c = —2 h ■ ([6. c]) with the segment [b,c] (cf. equation (3.13)), but because (3.12) does not hold 

b 

we will need to compensate for the discrepancy by adding the monopole term, p\ ,([6, c]) so that 


f(h) 


fi' <■]) = - 2 Jp* and A H) 


(M) = ^96r+p| r ([6,c]). 


(3.19) 


These considerations easily generalize to the entire grid boundary 7. Namely, for every pair of nodes (n \ , N2) 
such that ui G 7+ and n 2 G 7“ flN Ul , i.e., for every grid segment with the endpoints on different sides of T, 
we can introduce a discrete dipole with the moment (cf. equation (3.15)) 


Qn 1 


-2 h- 


/H ([»i,n 2 ]) 

In 1 


-2 


A^) 


Aft) 


h 



(3.20) 


and a discrete monopole with the intensity 

ft*([ni,»a]) = l [V/'.LM - ?<*>„»] • 

The reconstruction formulae that generalize (3.19) then become (cf. formulae (3.16)): 


v n 1 € 7 + : 4'* ) 
V n 2 € 7 _ : /* fe) 


E 2 h' qnin *' 


E 

m €'y + nN „ 2 


2/i Pll'2 ([^1 1 ^2]) 


(3.21) 


(3.22) 


We call the terms p U2 of the type (3.21) grid monopoles because they are due to “discontinuity” of the 
discrete derivative, i.e., discrepancy between the first difference quotients of and from the two sides 
of T, and they contribute toward the value of f\ h) only at one node of each respective pair (ni,tt 2 ), set 1 
(3.22). We must mention though that the latter decision, namely, to attribute the monopole sources from 
each appropriate grid segment only to the corresponding exterior endpoint (i.e., the endpoint n*> € 7“), is, of 
course, arbitrary. Instead, we could have attributed the necessary monopole correction (see formulae (3.18) 
and (3.19)) to the interior endpoint only (i.e., to n\ G 7 + ), or even distributed it in any proportion between 
the two endpoints ?ii and n 2 . Even though formally on the discrete level there is no clear advantage of 
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introducing monopoles one way or another, we will see (discussion after the proof of Lemma 3.3) that having 
the monopoles introduced the way we did it, i.e., only on 7 “, makes the representation of a given v w on 
N + in the form of a discrete single-layer potential most convenient. We should also note that if we assume 
boundedness of the first derivatives of v(x) and w(x) up to I\ and as such, boundedness of the corresponding 
first difference quotients (for small h), we conclude that both p n . 2 = 0(h ~ 1 ) and q llin . 2 - 0(h~ l ), whereas 
the monopole and dipole contributions to the actual grid density f-, hi have different orders with respect to 
the grid size 1 . 0{h~ l ) and 0(h~ 2 ), respectively, set* ( 3 . 22 ). 

Next, consider a special kind of namely, such that 


Vm, 67 : ui h J=v { n h J. (3.23) 

It is easy to see from (3.21) that — 0 , i.e., there will bo no dipole contribution to the overall density 

of the potential in this case. The only remaining contribution will be monopole, and one can rewrite (3.22) 
using (3.21) and (3.23) so that 


V Ui 6 7 + : /<*>■ = 0, 

1 Ini 


Vn,G 7 -: A h) \ n> =~l £ CW + < 1 

»3€{N-\7 + }nN M2 

E = Pn a - 


?ii + 


To obtain the the expression for f\ h} in (3.24), we have used the fact that Vm G 7 “: I = Q. 

«2 

which allowed us to equivalently replace the difference quotients of w {h} from ( 3 . 22 ) by the corresponding 
terms in (3.24). At the same time, this expression is, of course, an explicit formula for the application of the 
operator L {h) of (3.1) to the combined function u (h) of (3.4b). This formula obviously does not contain w { 
because is obtained from by truncation, as before. There is a minor difference, however, compared 
to the previously analyzed dipole case. Then, to obtain the expressions for dipole moments (3.15) we did 
need the boundary condition (3.14) that still involved w^ l + , even though f ~ h ! does not depend on in 
any event. Now. the boundary condition (3.23), which leads to monopoles p„ 2 , n 2 € 7~, see (3.24), does 
not contain the values either. As such, we can consider w^ h) only on N“\ 7 + throughout the entire 
argument when we summarize the construction of discrete monopoles in the following 

Definition 3.2. If the function is defined on N+ so that — 0 on M + , 

? IN+\7“ t 

and tc {li) = w^ l) is defined on N \7 + so that = 0 on M“\ 7 _ , and if boundary conditions (3-24) 

hold, then the corresponding f\ h) of (3,5) is called the density of the discrete single-layer potential , and 
expression (3.6) defines the discrete single-layer potential itself. 

Definition 3.2 is a natural finite- difference analogue to the notion of the continuous single-layer potential 
introduced in Section 2 . Indeed, a pure continuous single-layer potential is defined by the discontinuity 
in the derivative of the function [§^] r , but not in the function itself. Accordingly, formula ( 3 . 21 ) show's 
that the discrete monopole-type sources are driven by the discrepancies in the first difference quotients of 
the grid functions v {h) and w {h) on different sides of the interface. Moreover, boundary condition (3.23) 
that Definition 3.2 hinges upon is essentially a discrete Dirichlet boundary condition for the function w {h) . 
It says that there is no discontinuity in the function itself in this case, whereas the discontinuity in the 
derivative (more precisely, first difference quotients) is, generally speaking, allowed, and it drives the non- 
trivial potential density. Finally, the term “grid monopoles” is characteristic because the sources p U2 of 
(3.24) are concentrated only on 7 “, i.e., only on one layer of the grid boundary 7 . 
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Lemma 3.3. Any v (h) defined on N + such that L {h) v lh) = 0 on M + , where L {h) is the operator of (3.1) 
with c < 0, i.e. cither a Yukawa or Laplace operator, can be represented in the form of a discrete single-layer 
potential. 

Proof. What we essentially need to show is that for any v {h} that meets the requirements of the Lemma, 
there will always be a function w^ fi ^ on N _ \"y + that solves the homogeneous equation = 0 on 

M - \ 7 ~, and satisfies the far-field boundary condition (3.2) and boundary conditions (3.23) on 7 ". The 
problem of finding the appropriate function on N “\7 + can be classified as a discrete exterior problem 
of the Dirichlet type, because boundary conditions (3.23) prescribe the actual values of the discrete function 
w {h) on the boundary 7 “. In this sense, this problem is analogous to the continuous problem (2.14) that 
needs to be solved in order to represent a continuous i?(:c) on Q in the form of a single-layer potential. 

For the proof, we again employ a maximum principle argument, as we did when proving Lemma 3.2. In 
this case, though, it is completely straightforward. The corresponding homogeneous problem 

L {h) w {h) = 0 on MTVT, 

V 7*2 € 7 “ : = 0, 

obviously has as many equations as it has unknowns, and its only solution is trivial, because otherwise we 
would have had an interior-node positive maximum or negative minimum, which is not possible. □ 

Similarly to the previously analyzed double-layer case, for the Helmholtz equation the solvability of the 
foregoing discrete exterior Dirichlet problem for on N“\ 7 + needs to be assumed (this, again, means 
that the complementary domain is not resonant). In this cast 1 we can claim that a given v [h) on N + can be 
represented as a discrete single-layer potential. 

Let us also emphasize that because we defined the discrete single layer so that the monopole sources 
are concentrated only on 7 “, the corresponding potential (3.6) with the density } 1 l of (3.24) automatically 
reconstructs the function on the entire grid N + so that L {b h) {h) — 0 on all of the M + . As we will 
see in Section 4 , this makes the discrete single layer a borderline case of discrete volumetric controls that 
are derived for noise cancellation in the finite-difference framework. In contradistinction to that, previously 
(in the case of a double layer, as well as in all composition eases) the discrete surface potential would only 
reconstruct v {h) on N + \ 7 ~, and the values v ih) \~ will have to be complemented afterwards. 

To summarize, we shall sav that we have constructed discrete analogues of the continuous surface po- 
tentials for second-order operators. Similarly to the continuous case, the discrete potentials can be classified 
by the type of “discontinuity” of the grid function on the interface. The discrete potentials are driven by 
the densities that can be represented as combinations of grid monopoles and dipoles. Again, similarly to 
the continuous case, a given solution of the homogeneous finite-difference equation on the domain can be 
represented as a discrete surface potential in a variety of forms, including a pure single-layer potential, a pure 
double-layer potential, as well as combinations of the two, see formulae (3.22). For the Helmholtz equation, 
the aforementioned key result on representation of a given solution as either a single-layer potential or a 
double- layer potential require that the complementary domain be non-resonant. It is obvious, however, that 
resonances can always be avoided by changing the geometry of the complementary domain only. 

4. Discussion on Possible Applications. Assume that we have a given field v generated by some 
unknown sources located on Q\, and we want to eliminate it on II by active means, i.e., by adding the 
new source terms to the right-hand side of the overall governing equation so that their influence on Q will 
be exactly -v. As an important physical application in this context, one may think of active control of 
monochromatic sound described by the Helmholtz equation. In this case, the waves generated by sources on 
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n> should be treated as incoming with respect to the domain $1 (this field solves the homogeneous Helmholtz 
equation on 12), and the purpose of introducing the active control system is to cancel these waves out. A 
general solution to the noise control problem in terms of volumetric, rather than surface, control sources has 
been constructed in [9]. It is given by 


r(vol) __ r 

* control X/tC ^ i (4.1) 


where w is an auxiliary function that is supposed to satisfy only some relatively “loose” requirements, namely, 
u 6 L and w |j, — c| r , here v is the field that we want to control. The justification for formula 

(4.1) is iather straightforward. The field v to be controlled on f 1 can be annihilated by the surface potential 
with the density (anlr’^lr) 5 or “ ( Onlr* Hr) * w bich is the same. The latter can be equivalently 
rewritten on fl (using Green’s formula) as a volume potential with the density given by (4.1): 


V* € n : 


—v(x) = 


-w+ / QLwdy = ~ / QLwdy + / QLwdy 
Jn -h'i Jn 

- / SLudy = f 6f!:Z,dy. 

Jn, J d 


(4.2) 


One can also show that by considering the entire variety of auxiliary functions w that satisfy the aforemen- 
tioned requirements, we obtain all those and only those volumetric control sources that identically cancel 
the unwanted noise v on the domain 12. Moreover, the actual signal v that is needed to define the auxiliary 
function w (through the boundary values $£|,, and may, in fact, be “contaminated” in the sense that 
it may contain the component generated by the sources inside 12 (if any). This component is obviously a 
solution to the homogeneous equation on 12,, and shall be treated as an outgoing wave with respect to the 
domain 12. In the framework of [9] we rather consider the sound generated inside 1) as friendly, i.e., we do 
not want our control sources to have any effect on it. It turns out that in the presence of interior sound the 
formula (4.1) for controls does not change at all, even though the actual boundary values for w do change. 
Tin' reason is that the output of controls /<j™Jr 0 i gi ven bv (4.1) on 12 is the same as the surface potential with 
the density - (^| r ,w| r ) , see (4.2); and the latter, as we know from the analysis of Section 2, reconstructs 
on 12 only the incoming portion of the overall wave field, while being completely insensitive (on 12) to any 
outgoing component that may be present in the boundary data. 

In much the same spirit as the foregoing description of volumetric controls (see [9] for detail), Proposi- 
tion 2.1 of Section 2 describes the entire variety of surface controls that cancel a given field on the domain 
12. In this sense, the considerations of Section 2 complement those of [9], Obviously, the monopole and 
dipole surface sources described in Proposition 2.1 have to be taken with the minus sign for the purpose of 
eliminating a given field v(x), x € $2. Namely, 


0 

<H r ) - fa ([«’ - v)v6(D ) , (4.3) 

where w = w(x) is an auxiliary function, x G 12 u w e U, Lw = 0. Similarly to formula (4.1), the set 
of all appropriate surface controls, see (4.3), is parameterized by considering the entire variety of auxiliary 
functions w on 12 1 . Moreover, by its very construction the surface potential on 12 is insensitive to any field 
from the interior sources. In other words, if the actual outgoing waves are present they can be considered 
as an addition to the function-parameter u: . which has the meaning of an outgoing wave at any rati'. 


ir(surf) 

•'control 


dw dv 

On dn 
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Ill [10] we study different optimization formulations and strategies for the control sources (4.1) and (4.3). 
We emphasize that any control source from this class guarantees the identical cancellation of unwanted noise 
on the domain 0. Therefore, unlike in many approximate formulations (see. e.g., [8]), in which the extent 
of noise reduction is a key optimization criterion, in the current exact formulation it will not be a part 
of the function of merit. As such, we have rather looked at the criteria that are based primarily on the 
control sources themselves. For example, we have found that purely monopole surface controls from the 
class (4.3) obtained with w of (2.14) appear globally optimal (i.e., among both surface (4.3) and volumetric 
(4.1) controls) from the standpoint of minimizing the overall absolute strength, i.e., integral amplitude of 
the volume velocity (see [8] for the definition), of all sources. Mathematically, this translates to optimization 
in the sense of L\. At the same time, this surface monopole solution appears to radiate zero acoustic power 
(again, see [8] for the definition). There are other types of surface control sources that radiate no power. 
For example, the sources (4.3) that correspond to w = 0 radiate no power and also produce no reflection of 
the incoming wave to the domain fli. The latter circumstance, which is rather obvious in the framework of 
the analysis of this paper, has also been mentioned in [8]. On the other hand, purely dipole surface 1 controls 
(4.3) obtained with w of (2.15) are known to reflect all of the incoming wave back to the domain fti, see [8]. 
Among other combinations of type (4.3) one can, in fact, find those that would actually absorb the power of 
the incoming wave, i.e., be “beneficial” as far as energy. 4 Moreover, the extent of this power absorption can 
basically be made as large as one would want. This statement, which may look counterintuitive at a first 
glance, does not, of course, contradict the energy conservation. The explanation is that the field from the 
control sources (4.3) actually “loads” (see [8]) the original sources of sound located on fti, i.e., make them 
“push harder” against the medium and as such produce more energy, which is then partially consumed by 
the control system. 

In practice, of course, no control system can be designed using genuine continuous elements, as required 
by the previous analysis, e.g., continuous excitation (4.3) along the perimeter of the protected region. Actual 
technical devices, i.e., loudspeakers, that one will need in order to build such a system, have finite (often, 
small) size and can typically be assigned a point-wise location in space. As such, it is convenient to consider 
a discrete (finite-difference) formulation of the noise control problem from the very beginning. 

As shown in [10] (see also previous work [13 16]) a general solution for finite-difference volumetric control 
sources that cancel the unwanted noise on the grid subdomain is given by (cf. formula (4.1)) 


Ah, vol) 
J control 


-L (h 'w w I 


lM- 


(4.4) 


where js an auxiliary function that is defined on N” and is supposed to satisfy the selected far-field 
boundary condition (e.g., (3.2)) and the following boundary condition on 7: w {h) \„ ~ where v {fl] 

is the discrete field that we want to control. Other than that the function is arbitrary and as such, 
parameterizes the entire variety of appropriate volumetric control sources. The justification for formula 
(4.4) is based on the theory of generalized difference potentials, see [5,6]. Without elaborating on this issue 
here, we simply mention that once the auxiliary function w {h) meets the aforementioned requirements, the 
function G (h) f^t?o\ wil1 alwa T s coincide with ~v {h) on N+ , and as such, provide an ideal cancellation of 
the unwanted signal. Moreover, similarly to the continuous case, the actual boundary data v {h) \ that we 
need to define w w may actually contain a component due to the interior (i.e., outgoing) sound, if there are 1 
acoustic sources on M + . In this case, formula (4.4) does not change, and the resulting control sources will 

4 In the control theory such systems are said to possess the property of passivity, and they are known to be less likely to 
misbehave if the operational conditions appear to be out of range. 
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still be sensitive to the incoming component of the overall field only. This is dear because if we add to u^ h) 
any solution of the homogeneous equation defined on N“ , the control sources f ^ T ^ of (4.4) will not be 
affected. 

The discrete single-layer potential with the density (3.24) turns out to be a limiting, or borderline, case of 
the general solution (4.4). Indeed, let us define the function required by (4.4) as follows: | = , 

where is given, and = 0 on In other words, coincides with the given data on 

the grid boundary 7 and provides a solution to the homogeneous equation on the rest of the domain. This 
is obviously equivalent to solving the Dirichlet problem for w {h) on Vy+ with boundary conditions (3.23), 
and then complementing Having done that, we will obviously recover from (4.4) the same single 

layer as (3.24), only with the opposite sign for cancellation. 

In contradistinction to that, the discrete double-layer potential with the density (3.16), as well as any of 
the combination type discrete surface potentials with the density (3.22), is not a part of the general solution 
(4.4). It rather complements the general solution (4.4), in much the same way as surface* control sources 
(4.3) complemented the volumetric ones (4.1) in the continuous framework. The formula for discrete surface 
controls is obtained by taking (3.22) with the opposite sign 
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~2fi * qnyn - ~ 


(4.5) 


where the monopoles and dipoles are defined by (3.21) and (3.20), respectively, and dipoles are assumed not 
to be identically equal to zero, otherwise (4.5) will fall into the category (4.4). As one can easily conclude 
from the analysis of Section 3, discrete surface controls (4.5) are also insensitive to the outgoing sound; given 
the actual incoming field that we want to control, adding an outgoing component will simply alter the 
function-parameter t/4 h ) . 

By comparing formulae (4.4) and (4.5) one can conclude that from the standpoint of derivation the key 
distinction between the volumetric and surface control sources is rather formal: In (4.4) we do not allow 
any control sources on , whereas in (4.5) we do allow the control sources on the interior layer of the grid 
boundary 7 + and as such enlarge the overall family of the appropriate control sources. In practice, however, 
the benefit from introducing dipoles on 7 may be far more profound. It turns out that in many cases the 
engineering implementation of monopole sources is more difficult than that of dipole sources. As such, the 
discrete dipoles q ni1l2 may be introduced into a practical design as actual dipoles (oscillating double-sided 
membranes) located, e.g., at the centerpoints of the corresponding segments. 

Following the continuous case, we have also considered in [10] several optimization formulations for 
the discrete control sources. Clearly, the type of optimization problem that most easily lends itself to the 
numerical solution is minimization of the L > norm of the control sources (4.4). However, this optimization 
criterion lacks clear physical interpretation. As such, other criteria need to be employed. For example, 
similarly to the continuous case we have found that the* pure single layer on 7“ appears globally optimal 
(among all solutions (4.4)) in the sense of Li, which corresponds to minimizing the overall absolute strength 
of all sources. We have also found that the surface control sources (4.5) obtained for = 0 possess the 
nonreflecting property, i.e., they cancel out the unwanted noise on the protected region and at the same 
time do not alter the acoustic field on the complementary domain. This property on the discrete level 
follows immediately from the construction of the surface potential, in much the same way as it does on the 
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continuous level. For further detail on optimization of control sources for active cancellation of sound we 
refer the reader to our paper [10]. 

The last issue that yet remains to he addressed is the connection between the continuous and discrete 
formulations of the noise control problem. In other words, if we design a control system in the finite-difference 
framework, i.e., construct a discrete collection of noise-cancelling devices according to (4.5), (3.21), (3.20), 
what kind of performance shall we expect as far as eliminating the unwanted noise continuously throughout 
the protected region. First of all, to attempt a construction of the discrete control system, we need to 
require some standard properties from the discretization that we are using. Namely, the grid has to be 
sufficiently fine so that to well resolve the waves of the length A = 2 7r/k (recall, the operator L of (2.2) for 
c = k 2 becomes the Helmholtz operator so that equation Lu — f governs the propagation of weaves with the 
length A = 27 r/k). Moreover, the finite-difference scheme has to be consistent and stable (discretization (3.1) 
obviously satisfies the latter criteria). With these requirements met, the general theory of [5 7] basically 
says that the discrete potential v {h) defined on N + will approximate the continuous potential ??(#) defined 
on ft, provided that the discrete boundary data approximate the continuous boundary data (t—, e) in 
some special sense. Namely, once the continuous function and first-order normal derivative are known at 
the boundary T, normal derivatives of higher orders can be obtained via the differential equation itself, and 
the near-boundary values t;^ 1 can be calculated using Taylor’s expansion; the order of accuracy of the latter 
calculation with respect to h has to be at least as high as the order of accuracy of the interior scheme. In 
this case, the quality of the approximation, i.e., the rate of convergence of the discrete potential with respect 
to the grid size /?, will be the same as prescribed by the finite-difference scheme itself; this rate is 0(hr) for 
the particular example (3.1) that we have considered in Section 3. In other words, when designing an active 
control system following the finite-difference approach of Section 3, one can expect to have* the actual noise 
cancellation in the same approximate sense as the solution of the finite-difference equation approximates the 
corresponding solution of the original differential equation. 

Acknowledgment. The author is very thankful to Professor V. S. Ryaben kii for useful discussions on 
the subject of the paper. 
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